Nonlinear Regression (Archontoulis and Miguez) paper

https://femiguez.github.io/nlraa-docs/nlraa-Oddi-LFMC.html

https://onlinelibrary.wiley.com/doi/pdfdirect/10.1002/ece3.5543

# install.packages("emmeans")

library(nlme)
library(nlraa)
library(ggplot2)
library(emmeans) # emmeans::contrast(), emmeans::emmeans()
library(tidyverse)

Steps in fitting a nonlinear mixed model

Part I.

data(lfmc)
glimpse(lfmc)

Rows: 247
Columns: 6
$ leaf.type <fct> M. spinosum, M. spinosum, M. spinosum, S. bracteolactus, S. bracteolactus, S. bracteolactus, Gras…
$ time      <int> 1, 1, 1, 1, 1, 1, 1, 1, 13, 13, 13, 13, 13, 13, 13, 13, 13, 26, 26, 26, 26, 26, 26, 26, 26, 26, 3…
$ plot      <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ site      <fct> P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P…
$ lfmc      <dbl> 238.4, 269.8, 325.2, 326.7, 257.3, 279.3, 91.9, 95.0, 215.8, 216.7, 205.3, 208.5, 225.7, 246.8, 8…
$ group     <fct> SM plot 1, SM plot 1, SM plot 1, SS plot 1, SS plot 1, SS plot 1, GE plot 1, GE plot 1, SM plot 1…

lfmc

lfmc %>%
  count(leaf.type)

lfmc %>%
  count(time)

lfmc %>%
  count(plot)

lfmc %>%
  count(site)

lfmc %>%
  count(group)

lfmc %>%
  count(site, plot)

lfmc %>%
  ggplot(aes(time, lfmc)) +
  geom_point() +
  facet_wrap(~group)

lfmc.gd <- groupedData(lfmc ~ time | group, data = lfmc, order.groups = FALSE)
lfmc.gd

Grouped Data: lfmc ~ time | group
           leaf.type time plot site  lfmc     group
1        M. spinosum    1    1    P 238.4 SM plot 1
2        M. spinosum    1    1    P 269.8 SM plot 1
3        M. spinosum    1    1    P 325.2 SM plot 1
4   S. bracteolactus    1    1    P 326.7 SS plot 1
5   S. bracteolactus    1    1    P 257.3 SS plot 1
6   S. bracteolactus    1    1    P 279.3 SS plot 1
7            Grass E    1    1    P  91.9 GE plot 1
8            Grass E    1    1    P  95.0 GE plot 1
32       M. spinosum   13    1    P 215.8 SM plot 1
33       M. spinosum   13    1    P 216.7 SM plot 1
34       M. spinosum   13    1    P 205.3 SM plot 1
35  S. bracteolactus   13    1    P 208.5 SS plot 1
36  S. bracteolactus   13    1    P 225.7 SS plot 1
37  S. bracteolactus   13    1    P 246.8 SS plot 1
38           Grass E   13    1    P  81.3 GE plot 1
39           Grass E   13    1    P  82.4 GE plot 1
40           Grass E   13    1    P  78.1 GE plot 1
68       M. spinosum   26    1    P 194.4 SM plot 1
69       M. spinosum   26    1    P 166.0 SM plot 1
70       M. spinosum   26    1    P 212.0 SM plot 1
71  S. bracteolactus   26    1    P 232.2 SS plot 1
72  S. bracteolactus   26    1    P 190.1 SS plot 1
73  S. bracteolactus   26    1    P 241.6 SS plot 1
74           Grass E   26    1    P  53.9 GE plot 1
75           Grass E   26    1    P  71.9 GE plot 1
76           Grass E   26    1    P  63.1 GE plot 1
104      M. spinosum   36    1    P 178.0 SM plot 1
105      M. spinosum   36    1    P 165.9 SM plot 1
106      M. spinosum   36    1    P 144.0 SM plot 1
107 S. bracteolactus   36    1    P 144.1 SS plot 1
108 S. bracteolactus   36    1    P 142.7 SS plot 1
109 S. bracteolactus   36    1    P 170.0 SS plot 1
110          Grass E   36    1    P  45.7 GE plot 1
111          Grass E   36    1    P  37.3 GE plot 1
112          Grass E   36    1    P  33.3 GE plot 1
140      M. spinosum   57    1    P 103.5 SM plot 1
141      M. spinosum   57    1    P 115.2 SM plot 1
142      M. spinosum   57    1    P 115.2 SM plot 1
143 S. bracteolactus   57    1    P  76.7 SS plot 1
144 S. bracteolactus   57    1    P 110.5 SS plot 1
145 S. bracteolactus   57    1    P 148.5 SS plot 1
146          Grass E   57    1    P  23.0 GE plot 1
147          Grass E   57    1    P  23.3 GE plot 1
148          Grass E   57    1    P  25.6 GE plot 1
176      M. spinosum   68    1    P  65.8 SM plot 1
177      M. spinosum   68    1    P  86.9 SM plot 1
178      M. spinosum   68    1    P  79.5 SM plot 1
179 S. bracteolactus   68    1    P  70.3 SS plot 1
180 S. bracteolactus   68    1    P  67.5 SS plot 1
181 S. bracteolactus   68    1    P  87.6 SS plot 1
182          Grass E   68    1    P  14.0 GE plot 1
183          Grass E   68    1    P  24.3 GE plot 1
184          Grass E   68    1    P  17.1 GE plot 1
212      M. spinosum   80    1    P  67.0 SM plot 1
213      M. spinosum   80    1    P  70.5 SM plot 1
214      M. spinosum   80    1    P  57.9 SM plot 1
215 S. bracteolactus   80    1    P  70.6 SS plot 1
216 S. bracteolactus   80    1    P  55.4 SS plot 1
217 S. bracteolactus   80    1    P  86.0 SS plot 1
218          Grass E   80    1    P  11.4 GE plot 1
219          Grass E   80    1    P  16.2 GE plot 1
220          Grass E   80    1    P  12.8 GE plot 1
9        M. spinosum    1    2    P 251.7 SM plot 2
10       M. spinosum    1    2    P 238.7 SM plot 2
11       M. spinosum    1    2    P 252.6 SM plot 2
12  S. bracteolactus    1    2    P 257.6 SS plot 2
13  S. bracteolactus    1    2    P 265.7 SS plot 2
14  S. bracteolactus    1    2    P 273.9 SS plot 2
15           Grass E    1    2    P  60.6 GE plot 2
16           Grass E    1    2    P  53.7 GE plot 2
41       M. spinosum   13    2    P 228.0 SM plot 2
42       M. spinosum   13    2    P 231.4 SM plot 2
43       M. spinosum   13    2    P 226.6 SM plot 2
44  S. bracteolactus   13    2    P 222.7 SS plot 2
45  S. bracteolactus   13    2    P 243.9 SS plot 2
46  S. bracteolactus   13    2    P 218.6 SS plot 2
47           Grass E   13    2    P  50.5 GE plot 2
48           Grass E   13    2    P  64.3 GE plot 2
49           Grass E   13    2    P  64.3 GE plot 2
77       M. spinosum   26    2    P 185.8 SM plot 2
78       M. spinosum   26    2    P 169.7 SM plot 2
79       M. spinosum   26    2    P 154.8 SM plot 2
80  S. bracteolactus   26    2    P 207.2 SS plot 2
81  S. bracteolactus   26    2    P 200.0 SS plot 2
82  S. bracteolactus   26    2    P 204.9 SS plot 2
83           Grass E   26    2    P  50.5 GE plot 2
84           Grass E   26    2    P  43.2 GE plot 2
85           Grass E   26    2    P  52.7 GE plot 2
113      M. spinosum   36    2    P 115.9 SM plot 2
114      M. spinosum   36    2    P 124.4 SM plot 2
115      M. spinosum   36    2    P 130.7 SM plot 2
116 S. bracteolactus   36    2    P 214.9 SS plot 2
117 S. bracteolactus   36    2    P 159.2 SS plot 2
118 S. bracteolactus   36    2    P 135.9 SS plot 2
119          Grass E   36    2    P  32.1 GE plot 2
120          Grass E   36    2    P  34.0 GE plot 2
121          Grass E   36    2    P  39.2 GE plot 2
149      M. spinosum   57    2    P  96.0 SM plot 2
150      M. spinosum   57    2    P  76.7 SM plot 2
151      M. spinosum   57    2    P 101.2 SM plot 2
152 S. bracteolactus   57    2    P  80.6 SS plot 2
153 S. bracteolactus   57    2    P  83.7 SS plot 2
154 S. bracteolactus   57    2    P 108.0 SS plot 2
155          Grass E   57    2    P  20.4 GE plot 2
156          Grass E   57    2    P  23.7 GE plot 2
157          Grass E   57    2    P  21.7 GE plot 2
185      M. spinosum   68    2    P  62.2 SM plot 2
186      M. spinosum   68    2    P  70.6 SM plot 2
187      M. spinosum   68    2    P  76.6 SM plot 2
188 S. bracteolactus   68    2    P  72.8 SS plot 2
189 S. bracteolactus   68    2    P  85.2 SS plot 2
190 S. bracteolactus   68    2    P  75.4 SS plot 2
191          Grass E   68    2    P  15.4 GE plot 2
192          Grass E   68    2    P  11.8 GE plot 2
193          Grass E   68    2    P  15.3 GE plot 2
221      M. spinosum   80    2    P  38.1 SM plot 2
222      M. spinosum   80    2    P  61.2 SM plot 2
223      M. spinosum   80    2    P  53.2 SM plot 2
224 S. bracteolactus   80    2    P  72.2 SS plot 2
225 S. bracteolactus   80    2    P  64.2 SS plot 2
226 S. bracteolactus   80    2    P  61.4 SS plot 2
227          Grass E   80    2    P  13.9 GE plot 2
228          Grass E   80    2    P  11.2 GE plot 2
229          Grass E   80    2    P  13.3 GE plot 2
17       M. spinosum    1    3    P 224.3 SM plot 3
18       M. spinosum    1    3    P 207.3 SM plot 3
19  S. bracteolactus    1    3    P 240.3 SS plot 3
20  S. bracteolactus    1    3    P 215.8 SS plot 3
21           Grass E    1    3    P  75.3 GE plot 3
22           Grass E    1    3    P  65.7 GE plot 3
50       M. spinosum   13    3    P 240.8 SM plot 3
51       M. spinosum   13    3    P 178.4 SM plot 3
52       M. spinosum   13    3    P 231.2 SM plot 3
53  S. bracteolactus   13    3    P 218.6 SS plot 3
54  S. bracteolactus   13    3    P 244.1 SS plot 3
55  S. bracteolactus   13    3    P 204.4 SS plot 3
56           Grass E   13    3    P  69.3 GE plot 3
57           Grass E   13    3    P  77.2 GE plot 3
58           Grass E   13    3    P  62.5 GE plot 3
86       M. spinosum   26    3    P 154.7 SM plot 3
87       M. spinosum   26    3    P 221.0 SM plot 3
88       M. spinosum   26    3    P 231.7 SM plot 3
89  S. bracteolactus   26    3    P 171.4 SS plot 3
90  S. bracteolactus   26    3    P 209.4 SS plot 3
91  S. bracteolactus   26    3    P 180.5 SS plot 3
92           Grass E   26    3    P  54.7 GE plot 3
93           Grass E   26    3    P  45.9 GE plot 3
94           Grass E   26    3    P  41.1 GE plot 3
122      M. spinosum   36    3    P 164.7 SM plot 3
123      M. spinosum   36    3    P 142.9 SM plot 3
124      M. spinosum   36    3    P 161.4 SM plot 3
125 S. bracteolactus   36    3    P 152.0 SS plot 3
126 S. bracteolactus   36    3    P 114.6 SS plot 3
127 S. bracteolactus   36    3    P 183.8 SS plot 3
128          Grass E   36    3    P  24.4 GE plot 3
129          Grass E   36    3    P  28.0 GE plot 3
130          Grass E   36    3    P  26.4 GE plot 3
158      M. spinosum   57    3    P 129.1 SM plot 3
159      M. spinosum   57    3    P 119.9 SM plot 3
160      M. spinosum   57    3    P 148.1 SM plot 3
161 S. bracteolactus   57    3    P 110.6 SS plot 3
162 S. bracteolactus   57    3    P  54.1 SS plot 3
163 S. bracteolactus   57    3    P  95.1 SS plot 3
164          Grass E   57    3    P  31.3 GE plot 3
165          Grass E   57    3    P  22.5 GE plot 3
166          Grass E   57    3    P  27.4 GE plot 3
 [ reached 'max' / getOption("max.print") -- omitted 81 rows ]

lfmc.gd$plot <- with(lfmc.gd, factor(plot, levels = c("4", "5", "6", "1", "2", "3")))
lfmc.gd

Grouped Data: lfmc ~ time | group
           leaf.type time plot site  lfmc     group
1        M. spinosum    1    1    P 238.4 SM plot 1
2        M. spinosum    1    1    P 269.8 SM plot 1
3        M. spinosum    1    1    P 325.2 SM plot 1
4   S. bracteolactus    1    1    P 326.7 SS plot 1
5   S. bracteolactus    1    1    P 257.3 SS plot 1
6   S. bracteolactus    1    1    P 279.3 SS plot 1
7            Grass E    1    1    P  91.9 GE plot 1
8            Grass E    1    1    P  95.0 GE plot 1
32       M. spinosum   13    1    P 215.8 SM plot 1
33       M. spinosum   13    1    P 216.7 SM plot 1
34       M. spinosum   13    1    P 205.3 SM plot 1
35  S. bracteolactus   13    1    P 208.5 SS plot 1
36  S. bracteolactus   13    1    P 225.7 SS plot 1
37  S. bracteolactus   13    1    P 246.8 SS plot 1
38           Grass E   13    1    P  81.3 GE plot 1
39           Grass E   13    1    P  82.4 GE plot 1
40           Grass E   13    1    P  78.1 GE plot 1
68       M. spinosum   26    1    P 194.4 SM plot 1
69       M. spinosum   26    1    P 166.0 SM plot 1
70       M. spinosum   26    1    P 212.0 SM plot 1
71  S. bracteolactus   26    1    P 232.2 SS plot 1
72  S. bracteolactus   26    1    P 190.1 SS plot 1
73  S. bracteolactus   26    1    P 241.6 SS plot 1
74           Grass E   26    1    P  53.9 GE plot 1
75           Grass E   26    1    P  71.9 GE plot 1
76           Grass E   26    1    P  63.1 GE plot 1
104      M. spinosum   36    1    P 178.0 SM plot 1
105      M. spinosum   36    1    P 165.9 SM plot 1
106      M. spinosum   36    1    P 144.0 SM plot 1
107 S. bracteolactus   36    1    P 144.1 SS plot 1
108 S. bracteolactus   36    1    P 142.7 SS plot 1
109 S. bracteolactus   36    1    P 170.0 SS plot 1
110          Grass E   36    1    P  45.7 GE plot 1
111          Grass E   36    1    P  37.3 GE plot 1
112          Grass E   36    1    P  33.3 GE plot 1
140      M. spinosum   57    1    P 103.5 SM plot 1
141      M. spinosum   57    1    P 115.2 SM plot 1
142      M. spinosum   57    1    P 115.2 SM plot 1
143 S. bracteolactus   57    1    P  76.7 SS plot 1
144 S. bracteolactus   57    1    P 110.5 SS plot 1
145 S. bracteolactus   57    1    P 148.5 SS plot 1
146          Grass E   57    1    P  23.0 GE plot 1
147          Grass E   57    1    P  23.3 GE plot 1
148          Grass E   57    1    P  25.6 GE plot 1
176      M. spinosum   68    1    P  65.8 SM plot 1
177      M. spinosum   68    1    P  86.9 SM plot 1
178      M. spinosum   68    1    P  79.5 SM plot 1
179 S. bracteolactus   68    1    P  70.3 SS plot 1
180 S. bracteolactus   68    1    P  67.5 SS plot 1
181 S. bracteolactus   68    1    P  87.6 SS plot 1
182          Grass E   68    1    P  14.0 GE plot 1
183          Grass E   68    1    P  24.3 GE plot 1
184          Grass E   68    1    P  17.1 GE plot 1
212      M. spinosum   80    1    P  67.0 SM plot 1
213      M. spinosum   80    1    P  70.5 SM plot 1
214      M. spinosum   80    1    P  57.9 SM plot 1
215 S. bracteolactus   80    1    P  70.6 SS plot 1
216 S. bracteolactus   80    1    P  55.4 SS plot 1
217 S. bracteolactus   80    1    P  86.0 SS plot 1
218          Grass E   80    1    P  11.4 GE plot 1
219          Grass E   80    1    P  16.2 GE plot 1
220          Grass E   80    1    P  12.8 GE plot 1
9        M. spinosum    1    2    P 251.7 SM plot 2
10       M. spinosum    1    2    P 238.7 SM plot 2
11       M. spinosum    1    2    P 252.6 SM plot 2
12  S. bracteolactus    1    2    P 257.6 SS plot 2
13  S. bracteolactus    1    2    P 265.7 SS plot 2
14  S. bracteolactus    1    2    P 273.9 SS plot 2
15           Grass E    1    2    P  60.6 GE plot 2
16           Grass E    1    2    P  53.7 GE plot 2
41       M. spinosum   13    2    P 228.0 SM plot 2
42       M. spinosum   13    2    P 231.4 SM plot 2
43       M. spinosum   13    2    P 226.6 SM plot 2
44  S. bracteolactus   13    2    P 222.7 SS plot 2
45  S. bracteolactus   13    2    P 243.9 SS plot 2
46  S. bracteolactus   13    2    P 218.6 SS plot 2
47           Grass E   13    2    P  50.5 GE plot 2
48           Grass E   13    2    P  64.3 GE plot 2
49           Grass E   13    2    P  64.3 GE plot 2
77       M. spinosum   26    2    P 185.8 SM plot 2
78       M. spinosum   26    2    P 169.7 SM plot 2
79       M. spinosum   26    2    P 154.8 SM plot 2
80  S. bracteolactus   26    2    P 207.2 SS plot 2
81  S. bracteolactus   26    2    P 200.0 SS plot 2
82  S. bracteolactus   26    2    P 204.9 SS plot 2
83           Grass E   26    2    P  50.5 GE plot 2
84           Grass E   26    2    P  43.2 GE plot 2
85           Grass E   26    2    P  52.7 GE plot 2
113      M. spinosum   36    2    P 115.9 SM plot 2
114      M. spinosum   36    2    P 124.4 SM plot 2
115      M. spinosum   36    2    P 130.7 SM plot 2
116 S. bracteolactus   36    2    P 214.9 SS plot 2
117 S. bracteolactus   36    2    P 159.2 SS plot 2
118 S. bracteolactus   36    2    P 135.9 SS plot 2
119          Grass E   36    2    P  32.1 GE plot 2
120          Grass E   36    2    P  34.0 GE plot 2
121          Grass E   36    2    P  39.2 GE plot 2
149      M. spinosum   57    2    P  96.0 SM plot 2
150      M. spinosum   57    2    P  76.7 SM plot 2
151      M. spinosum   57    2    P 101.2 SM plot 2
152 S. bracteolactus   57    2    P  80.6 SS plot 2
153 S. bracteolactus   57    2    P  83.7 SS plot 2
154 S. bracteolactus   57    2    P 108.0 SS plot 2
155          Grass E   57    2    P  20.4 GE plot 2
156          Grass E   57    2    P  23.7 GE plot 2
157          Grass E   57    2    P  21.7 GE plot 2
185      M. spinosum   68    2    P  62.2 SM plot 2
186      M. spinosum   68    2    P  70.6 SM plot 2
187      M. spinosum   68    2    P  76.6 SM plot 2
188 S. bracteolactus   68    2    P  72.8 SS plot 2
189 S. bracteolactus   68    2    P  85.2 SS plot 2
190 S. bracteolactus   68    2    P  75.4 SS plot 2
191          Grass E   68    2    P  15.4 GE plot 2
192          Grass E   68    2    P  11.8 GE plot 2
193          Grass E   68    2    P  15.3 GE plot 2
221      M. spinosum   80    2    P  38.1 SM plot 2
222      M. spinosum   80    2    P  61.2 SM plot 2
223      M. spinosum   80    2    P  53.2 SM plot 2
224 S. bracteolactus   80    2    P  72.2 SS plot 2
225 S. bracteolactus   80    2    P  64.2 SS plot 2
226 S. bracteolactus   80    2    P  61.4 SS plot 2
227          Grass E   80    2    P  13.9 GE plot 2
228          Grass E   80    2    P  11.2 GE plot 2
229          Grass E   80    2    P  13.3 GE plot 2
17       M. spinosum    1    3    P 224.3 SM plot 3
18       M. spinosum    1    3    P 207.3 SM plot 3
19  S. bracteolactus    1    3    P 240.3 SS plot 3
20  S. bracteolactus    1    3    P 215.8 SS plot 3
21           Grass E    1    3    P  75.3 GE plot 3
22           Grass E    1    3    P  65.7 GE plot 3
50       M. spinosum   13    3    P 240.8 SM plot 3
51       M. spinosum   13    3    P 178.4 SM plot 3
52       M. spinosum   13    3    P 231.2 SM plot 3
53  S. bracteolactus   13    3    P 218.6 SS plot 3
54  S. bracteolactus   13    3    P 244.1 SS plot 3
55  S. bracteolactus   13    3    P 204.4 SS plot 3
56           Grass E   13    3    P  69.3 GE plot 3
57           Grass E   13    3    P  77.2 GE plot 3
58           Grass E   13    3    P  62.5 GE plot 3
86       M. spinosum   26    3    P 154.7 SM plot 3
87       M. spinosum   26    3    P 221.0 SM plot 3
88       M. spinosum   26    3    P 231.7 SM plot 3
89  S. bracteolactus   26    3    P 171.4 SS plot 3
90  S. bracteolactus   26    3    P 209.4 SS plot 3
91  S. bracteolactus   26    3    P 180.5 SS plot 3
92           Grass E   26    3    P  54.7 GE plot 3
93           Grass E   26    3    P  45.9 GE plot 3
94           Grass E   26    3    P  41.1 GE plot 3
122      M. spinosum   36    3    P 164.7 SM plot 3
123      M. spinosum   36    3    P 142.9 SM plot 3
124      M. spinosum   36    3    P 161.4 SM plot 3
125 S. bracteolactus   36    3    P 152.0 SS plot 3
126 S. bracteolactus   36    3    P 114.6 SS plot 3
127 S. bracteolactus   36    3    P 183.8 SS plot 3
128          Grass E   36    3    P  24.4 GE plot 3
129          Grass E   36    3    P  28.0 GE plot 3
130          Grass E   36    3    P  26.4 GE plot 3
158      M. spinosum   57    3    P 129.1 SM plot 3
159      M. spinosum   57    3    P 119.9 SM plot 3
160      M. spinosum   57    3    P 148.1 SM plot 3
161 S. bracteolactus   57    3    P 110.6 SS plot 3
162 S. bracteolactus   57    3    P  54.1 SS plot 3
163 S. bracteolactus   57    3    P  95.1 SS plot 3
164          Grass E   57    3    P  31.3 GE plot 3
165          Grass E   57    3    P  22.5 GE plot 3
166          Grass E   57    3    P  27.4 GE plot 3
 [ reached 'max' / getOption("max.print") -- omitted 81 rows ]

ggplot(lfmc, aes(time, lfmc)) +
  geom_point() +
  geom_smooth(method = "loess") +
  facet_wrap(~leaf.type) +
  xlab("Time (days)") +
  ylab("LMFC (%)")

`geom_smooth()` using formula 'y ~ x'

NONLINEAR MIXED-EFFECTS MODEL (M1)

nlsL <- nlsList(lfmc ~ SSdlf(time, A, w, m, s), data = lfmc.gd)

Warning: 3 errors caught in nls(model, data = data, control = controlvals).  The error messages and their frequencies are

step factor 0.000488281 reduced below 'minFactor' of 0.000976562 
                                                               1 
                     number of iterations exceeded maximum of 50 
                                                               2 

coef(nlsL)

plot(intervals(nlsL))

nl.re <- nlme(nlsL, control = nlmeControl(maxIter = 5000, msMaxIter = 1500))

Warning in (function (model, data = sys.frame(sys.parent()), fixed, random,  :
  Iteration 1, LME step: nlminb() did not converge (code = 1). PORT message: function evaluation limit reached without convergence (9)

nl.re.1 <- update(nl.re, random = pdDiag(A + w + m + s ~ 1))
fxf <- fixef(nl.re.1)
fxf

        A         w         m         s 
193.45267  21.45740  31.48270 -22.59711 

# nl.me <- update(nl.re.1, fixed = list(A + w + m + s ~ leaf.type),
#                 start = c(A = fxf[1],0,0,0,
#                           w = fxf[2],0,0,0,
#                           m = fxf[3],0,0,0,
#                           s = fxf[4],0,0,0))

nl.re.2 <- update(nl.re.1, random = pdDiag(A + w + m ~ 1))
nl.re.2

Nonlinear mixed-effects model fit by maximum likelihood
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 
  Log-likelihood: -1070.648
  Fixed: list(A ~ 1, w ~ 1, m ~ 1, s ~ 1) 
        A         w         m         s 
193.46656  21.43199  31.48851 -22.60602 

Random effects:
 Formula: list(A ~ 1, w ~ 1, m ~ 1)
 Level: group
 Structure: Diagonal
               A        w        m Residual
StdDev: 119.3792 10.94166 8.927276 15.26574

Number of Observations: 247
Number of Groups: 12 

fxf <- fixef(nl.re.2)
fxf

        A         w         m         s 
193.46656  21.43199  31.48851 -22.60602 

nl.me.1 <- update(
  nl.re.2,
  fixed = list(A + w + m + s ~ leaf.type),
  start = c(A = fxf[1], 0, 0, 0, w = fxf[2], 0, 0, 0, m = fxf[3], 0, 0, 0, s = fxf[4], 0, 0, 0)
)
nl.me.1

Nonlinear mixed-effects model fit by maximum likelihood
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 
  Log-likelihood: -1034.375
  Fixed: list(A + w + m + s ~ leaf.type) 
              A.(Intercept)          A.leaf.typeGrass E      A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus 
                  23.135414                   53.294908                  296.620596                  263.313118 
              w.(Intercept)          w.leaf.typeGrass E      w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus 
                  50.291577                  -34.313924                  -26.732287                    2.135665 
              m.(Intercept)          m.leaf.typeGrass E      m.leaf.typeM. spinosum m.leaf.typeS. bracteolactus 
                  51.533983                  -21.920404                  -20.647403                  -18.170359 
              s.(Intercept)          s.leaf.typeGrass E      s.leaf.typeM. spinosum s.leaf.typeS. bracteolactus 
                  17.605794                  -25.926031                  -43.529450                  -33.666226 

Random effects:
 Formula: list(A ~ 1, w ~ 1, m ~ 1)
 Level: group
 Structure: Diagonal
        A.(Intercept) w.(Intercept) m.(Intercept) Residual
StdDev:      8.112813   0.001236897      2.468334 15.38607

Number of Observations: 247
Number of Groups: 12 

AIC(nl.re.2, nl.me.1)

IC_tab(nl.re.2, nl.me.1)

The small StdDev of “w” suggests to remove the random-effect of plot on this parameter.

nl.me.2 <- update(nl.me.1, random = pdDiag(A + m ~ 1))
IC_tab(nl.me.1, nl.me.2)

in addition, if the random-effects on “m” is also removed.

nl.me.3 <- update(nl.me.2, random = pdDiag(A ~ 1))
nl.me.3

Nonlinear mixed-effects model fit by maximum likelihood
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 
  Log-likelihood: -1035.19
  Fixed: list(A + w + m + s ~ leaf.type) 
              A.(Intercept)          A.leaf.typeGrass E      A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus 
                  23.071500                   53.950821                  408.256907                  258.698322 
              w.(Intercept)          w.leaf.typeGrass E      w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus 
                  50.311106                  -34.625804                  -79.299174                    3.717378 
              m.(Intercept)          m.leaf.typeGrass E      m.leaf.typeM. spinosum m.leaf.typeS. bracteolactus 
                  51.608221                  -21.975396                  -32.109785                  -17.665025 
              s.(Intercept)          s.leaf.typeGrass E      s.leaf.typeM. spinosum s.leaf.typeS. bracteolactus 
                  17.683236                  -26.518320                  -60.147105                  -33.113726 

Random effects:
 Formula: A ~ 1 | group
        A.(Intercept) Residual
StdDev:       10.6035 15.53457

Number of Observations: 247
Number of Groups: 12 

the model fit is similar, supporting the simpler model.

IC_tab(nl.me.2, nl.me.3)

plot(nl.me.3)

hist(residuals(nl.me.3, type = "normalized"), breaks = 20)

rse <- residuals(nl.me.3) / sigma(nl.me.3)
par(pty = "s")
qqnorm(rse)
qqline(rse)

Heterocedasticity: variance modeling —-

nl.me.3.vm <- update(nl.me.3, weights = varIdent(form = ~ 1 | leaf.type))
nl.me.3.vm

Nonlinear mixed-effects model fit by maximum likelihood
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 
  Log-likelihood: -965.5628
  Fixed: list(A + w + m + s ~ leaf.type) 
              A.(Intercept)          A.leaf.typeGrass E      A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus 
                 50.1522851                  28.9553861                 346.5816334                 231.1569667 
              w.(Intercept)          w.leaf.typeGrass E      w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus 
                 24.6693607                  -9.4837876                 -35.8327942                  29.7533014 
              m.(Intercept)          m.leaf.typeGrass E      m.leaf.typeM. spinosum m.leaf.typeS. bracteolactus 
                 48.9613646                 -19.9780200                 -26.2920794                 -14.9909275 
              s.(Intercept)          s.leaf.typeGrass E      s.leaf.typeM. spinosum s.leaf.typeS. bracteolactus 
                -16.2294127                   6.3520529                 -21.0664455                   0.9389323 

Random effects:
 Formula: A ~ 1 | group
        A.(Intercept) Residual
StdDev:      9.412263 6.824811

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | leaf.type 
 Parameter estimates:
         Grass W          Grass E      M. spinosum S. bracteolactus 
       1.0000000        0.8547835        3.1456077        3.0262116 
Number of Observations: 247
Number of Groups: 12 

IC_tab(nl.me.3, nl.me.3.vm)

여기서 가중잔차는 weight를 이용해서 계산해야 제대로 plot이 됨

res <-
  tibble(
    fitted = fitted(nl.me.3.vm),
    resid = resid(nl.me.3.vm), # 일반 잔차
    nresid = resid(nl.me.3.vm, type = "normalized"), # 표준 잔차
    weight = diag(var_cov(nl.me.3.vm)), # 가중치
    wresid = sqrt(1 / weight) * resid # 가중잔차 수동계산
  ) %>%
  bind_cols(lfmc.gd)
res

res %>%
  ggplot(aes(fitted, nresid)) +
  geom_point(aes(color = leaf.type))

res %>%
  ggplot(aes(fitted, wresid)) +
  geom_point(aes(color = leaf.type))

plot(nl.me.3.vm)

hist(residuals(nl.me.3.vm, type = "normalized"), main = "", xlab = "Residuals")

the residuals look fine now.

qqnorm(resid(nl.me.3.vm, type = "normalized"))
qqline(resid(nl.me.3.vm, type = "normalized"))

Temporal correlation check. This plot suggests there is no need for modeling the correlation of the residuals.

plot(ACF(nl.me.3.vm), alpha = 0.01)

Evaluation of Fixed Effects

the model is fitted by Maximum likelihood

nl.me.ml <- update(nl.me.3.vm, method = "ML")
nl.me.ml

Nonlinear mixed-effects model fit by maximum likelihood
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 
  Log-likelihood: -965.5628
  Fixed: list(A + w + m + s ~ leaf.type) 
              A.(Intercept)          A.leaf.typeGrass E      A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus 
                 50.1522851                  28.9553861                 346.5816334                 231.1569667 
              w.(Intercept)          w.leaf.typeGrass E      w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus 
                 24.6693607                  -9.4837876                 -35.8327942                  29.7533014 
              m.(Intercept)          m.leaf.typeGrass E      m.leaf.typeM. spinosum m.leaf.typeS. bracteolactus 
                 48.9613646                 -19.9780200                 -26.2920794                 -14.9909275 
              s.(Intercept)          s.leaf.typeGrass E      s.leaf.typeM. spinosum s.leaf.typeS. bracteolactus 
                -16.2294127                   6.3520529                 -21.0664455                   0.9389323 

Random effects:
 Formula: A ~ 1 | group
        A.(Intercept) Residual
StdDev:      9.412263 6.824811

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | leaf.type 
 Parameter estimates:
         Grass W          Grass E      M. spinosum S. bracteolactus 
       1.0000000        0.8547835        3.1456077        3.0262116 
Number of Observations: 247
Number of Groups: 12 

the model is reduced removing the fixed-effect of “leaf type” on “s”

nl.me.ml.1 <- update(nl.me.ml,
  fixed = list(A + w + m ~ leaf.type, s ~ 1),
  start = c(A = fxf[1], 0, 0, 0, w = fxf[2], 0, 0, 0, m = fxf[3], 0, 0, 0, s = fxf[4])
)
nl.me.ml.1

Nonlinear mixed-effects model fit by maximum likelihood
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 
  Log-likelihood: -969.2617
  Fixed: list(A + w + m ~ leaf.type, s ~ 1) 
              A.(Intercept)          A.leaf.typeGrass E      A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus 
                   49.94562                    39.60979                   219.45652                   232.11909 
              w.(Intercept)          w.leaf.typeGrass E      w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus 
                   25.24501                   -14.34293                    35.36122                    28.70389 
              m.(Intercept)          m.leaf.typeGrass E      m.leaf.typeM. spinosum m.leaf.typeS. bracteolactus 
                   48.30390                   -21.59477                   -15.57765                   -14.37454 
                          s 
                  -15.44291 

Random effects:
 Formula: A ~ 1 | group
        A.(Intercept) Residual
StdDev:      9.612472 6.825604

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | leaf.type 
 Parameter estimates:
         Grass W          Grass E      M. spinosum S. bracteolactus 
       1.0000000        0.8714808        3.2750696        3.0238816 
Number of Observations: 247
Number of Groups: 12 

the AIC comparison suggests that the more complex model is not justified (fits are similar)

IC_tab(nl.me.ml, nl.me.ml.1)

the model is reduced removing the fixed-effect of “leaf type” on “m”

nl.me.ml.2 <- update(nl.me.ml,
  fixed = list(A + w ~ leaf.type, m + s ~ 1),
  start = c(A = fxf[1], 0, 0, 0, w = fxf[2], 0, 0, 0, m = fxf[3], s = fxf[4])
)
nl.me.ml.2

Nonlinear mixed-effects model fit by maximum likelihood
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 
  Log-likelihood: -972.7465
  Fixed: list(A + w ~ leaf.type, m + s ~ 1) 
              A.(Intercept)          A.leaf.typeGrass E      A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus 
                   54.25730                    30.93396                   223.27914                   240.31542 
              w.(Intercept)          w.leaf.typeGrass E      w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus 
                   29.05803                   -20.68663                    31.68646                    26.99003 
                          m                           s 
                   30.91655                   -16.15435 

Random effects:
 Formula: A ~ 1 | group
        A.(Intercept) Residual
StdDev:      9.216113 7.028345

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | leaf.type 
 Parameter estimates:
         Grass W          Grass E      M. spinosum S. bracteolactus 
       1.0000000        0.8820414        3.1693973        2.9643303 
Number of Observations: 247
Number of Groups: 12 

the model is reduced removing the fixed-effect of “leaf type” on “W”.

nl.me.ml.3 <- update(nl.me.ml,
  fixed = list(A ~ leaf.type, w + m + s ~ 1),
  start = c(A = fxf[1], 0, 0, 0, w = fxf[2], m = fxf[3], s = fxf[4])
)
nl.me.ml.3

Nonlinear mixed-effects model fit by maximum likelihood
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 
  Log-likelihood: -1017.91
  Fixed: list(A ~ leaf.type, w + m + s ~ 1) 
              A.(Intercept)          A.leaf.typeGrass E      A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus 
                   63.92200                    14.27252                   251.99435                   265.42764 
                          w                           m                           s 
                   17.33019                    32.65673                   -25.16298 

Random effects:
 Formula: A ~ 1 | group
        A.(Intercept) Residual
StdDev:      8.138947 7.795662

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | leaf.type 
 Parameter estimates:
         Grass W          Grass E      M. spinosum S. bracteolactus 
        1.000000         1.584887         2.813985         2.731415 
Number of Observations: 247
Number of Groups: 12 

nl.me.ml.4 <- update(nl.me.ml,
  fixed = list(w ~ leaf.type, A + m + s ~ 1),
  start = c(A = fxf[1], w = fxf[2], 0, 0, 0, m = fxf[3], s = fxf[4])
)
nl.me.ml.4

Nonlinear mixed-effects model fit by maximum likelihood
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 
  Log-likelihood: -999.8224
  Fixed: list(w ~ leaf.type, A + m + s ~ 1) 
              w.(Intercept)          w.leaf.typeGrass E      w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus 
                   29.16029                   -20.38028                    33.19430                    28.72481 
                          A                           m                           s 
                  175.78131                    31.08598                   -15.70297 

Random effects:
 Formula: A ~ 1 | group
               A Residual
StdDev: 108.0655 7.042308

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | leaf.type 
 Parameter estimates:
         Grass W          Grass E      M. spinosum S. bracteolactus 
       1.0000000        0.8770721        3.1889180        2.9304615 
Number of Observations: 247
Number of Groups: 12 

IC_tab(nl.me.ml.1, nl.me.ml.2, nl.me.ml.3, nl.me.ml.4)

Final Model

This model assumes:

• “A” and “w” vary with leaf type (fixed-effects).
• “A” varies with plot (random-effects).
• “m” and “s” are unique for all the plots and leaf types.
• Residual variance depends on “leaf type”
• There is no temporal dependence
• Residuals following a normal distribution

M1 <- update(nl.me.ml.2, method = "REML")
summary(M1)

Nonlinear mixed-effects model fit by REML
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 
       AIC      BIC   logLik
  1931.668 1983.689 -950.834

Random effects:
 Formula: A ~ 1 | group
        A.(Intercept) Residual
StdDev:      9.408521 7.162899

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | leaf.type 
 Parameter estimates:
         Grass W          Grass E      M. spinosum S. bracteolactus 
       1.0000000        0.8851764        3.1796102        2.9647393 
Fixed effects:  list(A + w ~ leaf.type, m + s ~ 1) 
                                Value Std.Error  DF   t-value p-value
A.(Intercept)                54.25350  5.988181 226  9.060097   0e+00
A.leaf.typeGrass E           30.92293  8.835210 226  3.499966   6e-04
A.leaf.typeM. spinosum      223.24504 15.915557 226 14.026844   0e+00
A.leaf.typeS. bracteolactus 240.27654 16.857356 226 14.253513   0e+00
w.(Intercept)                29.05581  1.712960 226 16.962334   0e+00
w.leaf.typeGrass E          -20.68938  2.592469 226 -7.980569   0e+00
w.leaf.typeM. spinosum       31.67297  7.753176 226  4.085161   1e-04
w.leaf.typeS. bracteolactus  26.97508  8.071111 226  3.342177   1e-03
m                            30.92881  2.065196 226 14.976213   0e+00
s                           -16.15406  2.365392 226 -6.829339   0e+00
 Correlation: 
                            A.(In) A.l.GE A..M.s A..S.b w.(In) w.l.GE w..M.s w..S.b m     
A.leaf.typeGrass E          -0.527                                                        
A.leaf.typeM. spinosum      -0.146  0.535                                                 
A.leaf.typeS. bracteolactus -0.115  0.537  0.746                                          
w.(Intercept)               -0.217 -0.034 -0.200 -0.216                                   
w.leaf.typeGrass E          -0.035 -0.283 -0.382 -0.399 -0.258                            
w.leaf.typeM. spinosum      -0.118 -0.227 -0.547 -0.454  0.157  0.573                     
w.leaf.typeS. bracteolactus -0.129 -0.241 -0.462 -0.575  0.188  0.601  0.640              
m                           -0.217 -0.324 -0.633 -0.666  0.092  0.145  0.161  0.172       
s                           -0.246 -0.354 -0.710 -0.748  0.416  0.575  0.702  0.752  0.544

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.27999661 -0.59657156 -0.02042043  0.57273326  3.12018001 

Number of Observations: 247
Number of Groups: 12 

fixef(M1)

              A.(Intercept)          A.leaf.typeGrass E      A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus 
                   54.25350                    30.92293                   223.24504                   240.27654 
              w.(Intercept)          w.leaf.typeGrass E      w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus 
                   29.05581                   -20.68938                    31.67297                    26.97508 
                          m                           s 
                   30.92881                   -16.15406 

ranef(M1)

          A.(Intercept)
GW plot 4     -1.915443
GW plot 5     -1.923198
GW plot 6      3.838641
GE plot 1     15.652459
SM plot 1      6.967303
SS plot 1      9.166925
GE plot 2    -11.066197
SM plot 2     -5.358910
SS plot 2      2.442781
GE plot 3     -4.586263
SM plot 3     -1.608393
SS plot 3    -11.609707

intervals(M1)

Approximate 95% confidence intervals

 Fixed effects:
                                lower      est.     upper
A.(Intercept)                42.45370  54.25350  66.05331
A.leaf.typeGrass E           13.51301  30.92293  48.33286
A.leaf.typeM. spinosum      191.88318 223.24504 254.60691
A.leaf.typeS. bracteolactus 207.05884 240.27654 273.49423
w.(Intercept)                25.68039  29.05581  32.43122
w.leaf.typeGrass E          -25.79788 -20.68938 -15.58088
w.leaf.typeM. spinosum       16.39521  31.67297  46.95073
w.leaf.typeS. bracteolactus  11.07083  26.97508  42.87934
m                            26.85931  30.92881  34.99831
s                           -20.81510 -16.15406 -11.49302

 Random Effects:
  Level: group 
                   lower     est.    upper
sd(A.(Intercept)) 5.0352 9.408521 17.58029

 Variance function:
                    lower      est.    upper
Grass E          0.675418 0.8851764 1.160077
M. spinosum      2.451780 3.1796102 4.123502
S. bracteolactus 2.285695 2.9647393 3.845517

 Within-group standard error:
   lower     est.    upper 
5.966409 7.162899 8.599330 

Part II

# 2.2.1 - Overall predictions (Table 3) ----
Agw <- fixef(M1)[1] # "A" for grasses in the W site (GW)
Age <- fixef(M1)[1] + fixef(M1)[2] # "A" for grasses in the E site (GE)
Asm <- fixef(M1)[1] + fixef(M1)[3] # "A" for Mullinum spinosum (SM)
Ass <- fixef(M1)[1] + fixef(M1)[4] # "A" for Senecio filaginoides (SS)
Wgw <- fixef(M1)[5] # "w" for grasses in the W site (GW)
Wge <- fixef(M1)[5] + fixef(M1)[6] # "w" for grasses in the E site (GE)
Wsm <- fixef(M1)[5] + fixef(M1)[7] # "w" for Mullinum spinosum (SM)
Wss <- fixef(M1)[5] + fixef(M1)[8] # "w" for Senecio filaginoides (SS)
M <- fixef(M1)[9] # "m" for all the leaf types
S <- fixef(M1)[10] # "s" for all the leaf types

GW <- subset(lfmc, leaf.type == "Grass W")
GE <- subset(lfmc, leaf.type == "Grass E")
SM <- subset(lfmc, leaf.type == "M. spinosum")
SS <- subset(lfmc, leaf.type == "S. bracteolactus")

par(mfrow = c(2, 2), cex = 0.75)
# Grasses W
plot(GW$time, GW$lfmc, xlab = "", ylab = "LFMC (%)", ylim = c(0, 100), main = "Grasses W site", col = "darkgreen")
curve(((Agw - Wgw) / (1 + exp((M - x) / S)) + Wgw), lwd = 2, lty = 1, add = T, col = "darkgreen")

# Grasses E
plot(GE$time, GE$lfmc, xlab = "", ylab = "LFMC (%)", ylim = c(0, 100), main = "Grasses E site", col = "green")
curve(((Age - Wge) / (1 + exp((M - x) / S)) + Wge), lwd = 2, lty = 1, add = T, col = "green")

# M. Spinosum
plot(SM$time, SM$lfmc, xlab = "Time (days)", ylab = "LFMC (%)", ylim = c(0, 350), main = "M. spinosum", font.main = 4, col = "darkorange")
curve(((Asm - Wsm) / (1 + exp((M - x) / S)) + Wsm), lwd = 2, lty = 1, add = T, col = "darkorange")

# S. filaginoides
plot(SS$time, SS$lfmc, xlab = "Time (days)", ylab = "LFMC (%)", ylim = c(0, 350), main = "S. bracteolactus", font.main = 4, col = "darkred")
curve(((Ass - Wss) / (1 + exp((M - x) / S)) + Wss), lwd = 2, lty = 1, add = T, col = "darkred")

ranef(M1)

          A.(Intercept)
GW plot 4     -1.915443
GW plot 5     -1.923198
GW plot 6      3.838641
GE plot 1     15.652459
SM plot 1      6.967303
SS plot 1      9.166925
GE plot 2    -11.066197
SM plot 2     -5.358910
SS plot 2      2.442781
GE plot 3     -4.586263
SM plot 3     -1.608393
SS plot 3    -11.609707

# 2.2.2 - Predictions at plot level (Table 3) ----
Agw.p4 <- fixef(M1)[1] + ranef(M1)[1, 1] # "A" for grasses in the plot 4 (GW)
Agw.p5 <- fixef(M1)[1] + ranef(M1)[2, 1] # "A" for grasses in the plot 5 (GW)
Agw.p6 <- fixef(M1)[1] + ranef(M1)[3, 1] # "A" for grasses in the plot 6 (GW)
Age.p1 <- fixef(M1)[1] + fixef(M1)[2] + ranef(M1)[4, 1] # "A" for grasses in the plot 1 (GW)
Age.p2 <- fixef(M1)[1] + fixef(M1)[2] + ranef(M1)[7, 1] # "A" for grasses in the plot 2 (GW)
Age.p3 <- fixef(M1)[1] + fixef(M1)[2] + ranef(M1)[10, 1] # "A" for grasses in the plot 3 (GW)
Asm.p1 <- fixef(M1)[1] + fixef(M1)[3] + ranef(M1)[5, 1] # "A" for M. spinosum in the plot 1 (SM)
Asm.p2 <- fixef(M1)[1] + fixef(M1)[3] + ranef(M1)[8, 1] # "A" for M. spinosum in the plot 2 (SM)
Asm.p3 <- fixef(M1)[1] + fixef(M1)[3] + ranef(M1)[11, 1] # "A" for M. spinosum in the plot 3 (SM)
Ass.p1 <- fixef(M1)[1] + fixef(M1)[4] + ranef(M1)[6, 1] # "A" for S. filaginoides in the plot 1 (SM)
Ass.p2 <- fixef(M1)[1] + fixef(M1)[4] + ranef(M1)[9, 1] # "A" for S. filaginoides in the plot 2 (SM)
Ass.p3 <- fixef(M1)[1] + fixef(M1)[4] + ranef(M1)[12, 1] # "A" for S. filaginoides in the plot 3 (SM)

par(mfrow = c(2, 2), cex = 0.75) # (Figure 4)
# Grasses W site
plot(GW$time, GW$lfmc, xlab = "", ylab = "LFMC (%)", ylim = c(0, 100), main = "Grasses W site", col = "darkgreen")
curve(((Agw.p4 - Wgw) / (1 + exp((M - x) / S)) + Wgw), lwd = 1, lty = 2, add = T, col = "darkgreen") # plot 4

curve(((Agw.p5 - Wgw) / (1 + exp((M - x) / S)) + Wgw), lwd = 1, lty = 2, add = T, col = "darkgreen") # plot 5
curve(((Agw.p6 - Wgw) / (1 + exp((M - x) / S)) + Wgw), lwd = 1, lty = 2, add = T, col = "darkgreen") # plot 6

# Grasses E site
plot(GE$time, GE$lfmc, xlab = "", ylab = "LFMC (%)", ylim = c(0, 100), main = "Grasses E site", col = "green")
curve(((Age.p1 - Wge) / (1 + exp((M - x) / S)) + Wge), lwd = 1, lty = 2, add = T, col = "green") # plot 1

curve(((Age.p2 - Wge) / (1 + exp((M - x) / S)) + Wge), lwd = 1, lty = 2, add = T, col = "green") # plot 2
curve(((Age.p3 - Wge) / (1 + exp((M - x) / S)) + Wge), lwd = 1, lty = 2, add = T, col = "green") # plot 3

# M. Spinosum (E site)
plot(SM$time, SM$lfmc, xlab = "Time (days)", ylab = "LFMC (%)", ylim = c(0, 350), main = "M. spinosum", font.main = 4, col = "darkorange")
curve(((Asm.p1 - Wsm) / (1 + exp((M - x) / S)) + Wsm), lwd = 1, lty = 2, add = T, col = "darkorange") # plot 1

curve(((Asm.p2 - Wsm) / (1 + exp((M - x) / S)) + Wsm), lwd = 1, lty = 2, add = T, col = "darkorange") # plot 2
curve(((Asm.p3 - Wsm) / (1 + exp((M - x) / S)) + Wsm), lwd = 1, lty = 2, add = T, col = "darkorange") # plot 3

# S. filaginoides (E site)
plot(SS$time, SS$lfmc, xlab = "Time (days)", ylab = "LFMC (%)", ylim = c(0, 350), main = "S. bracteolactus", font.main = 4, col = "darkred")
curve(((Ass.p1 - Wss) / (1 + exp((M - x) / S)) + Wss), lwd = 1, lty = 2, add = T, col = "darkred") # plot 1

curve(((Ass.p2 - Wss) / (1 + exp((M - x) / S)) + Wss), lwd = 1, lty = 2, add = T, col = "darkred") # plot 2
curve(((Ass.p3 - Wss) / (1 + exp((M - x) / S)) + Wss), lwd = 1, lty = 2, add = T, col = "darkred") # plot 3

# 2.2.3 - Derivatives (drying speed) ----
lfmc.GW <- expression((Agw - Wgw) / (1 + exp((M - x) / S)) + Wgw) # LFMC dynamics for grasses in the W site
ds.GW <- deriv(lfmc.GW, "x", function.arg = TRUE) # Drying speed for grasses in the W site

lfmc.GE <- expression((Age - Wge) / (1 + exp((M - x) / S)) + Wge) # LFMC dynamics for grasses in the E site
ds.GE <- deriv(lfmc.GE, "x", function.arg = TRUE) # Drying speed for grasses in the E site

lfmc.SM <- expression((Asm - Wsm) / (1 + exp((M - x) / S)) + Wsm) # LFMC dynamics for M. spinosum
ds.SM <- deriv(lfmc.SM, "x", function.arg = TRUE) # Drying speed for M. spinosum

lfmc.SS <- expression((Ass - Wss) / (1 + exp((M - x) / S)) + Wss) # LFMC dynamics for S. filaginoides
ds.SS <- deriv(lfmc.SS, "x", function.arg = TRUE) # # Drying speed for S. filaginoides

xvec <- seq(0, 100, length = 1000)
y.ds.GW <- ds.GW(xvec)
y.ds.GE <- ds.GE(xvec)
y.ds.SM <- ds.SM(xvec)
y.ds.SS <- ds.SS(xvec)

par(mfrow = c(2, 2), cex = 0.75) # (Figure 4)
plot(xvec, y.ds.GW,
  xlim = c(0, 80), ylim = c(0, 4),
  ylab = "Drying speed", xlab = "Time (days)", type = "n", main = "Grasses in the W site"
)
lines(xvec, -1 * (attr(y.ds.GW, "grad")), lwd = 2, lty = 1, col = "darkgreen")

plot(xvec, y.ds.GE,
  xlim = c(0, 80), ylim = c(0, 4),
  ylab = "Drying speed", xlab = "Time (days)", type = "n", main = "Grasses in the E site"
)
lines(xvec, -1 * (attr(y.ds.GE, "grad")), lwd = 2, lty = 1, col = "green")

plot(xvec, y.ds.SM,
  xlim = c(0, 80), ylim = c(0, 4),
  ylab = "Drying speed", xlab = "Time (days)", type = "n", main = "M. spinosum", font.main = 4
)
lines(xvec, -1 * (attr(y.ds.SM, "grad")), lwd = 2, lty = 1, col = "darkorange")

plot(xvec, y.ds.SS,
  xlim = c(0, 80), ylim = c(0, 4),
  ylab = "Drying speed", xlab = "Time (days)", type = "n", main = "S. filaginoides", font.main = 4
)
lines(xvec, -1 * (attr(y.ds.SS, "grad")), lwd = 2, lty = 1, col = "darkred")

# 2.3 ---- Testing differences among leaf types -------------------------------------------

# 2.3.1 - Parameter "A" ----
contrast(emmeans(M1, ~leaf.type, param = "A"), "pairwise")

 contrast                       estimate    SE  df t.ratio p.value
 Grass W - Grass E                 -30.9  8.84 226  -3.500  0.0031
 Grass W - M. spinosum            -223.2 15.92 226 -14.027  <.0001
 Grass W - S. bracteolactus       -240.3 16.86 226 -14.254  <.0001
 Grass E - M. spinosum            -192.3 13.45 226 -14.294  <.0001
 Grass E - S. bracteolactus       -209.4 14.22 226 -14.718  <.0001
 M. spinosum - S. bracteolactus    -17.0 11.71 226  -1.454  0.4671

P value adjustment: tukey method for comparing a family of 4 estimates 

• Maximum LFMC in grasses from the W site is lower than grasses from the E site (Grass W - Grass E)
• Maximum LFMC in grasses is lower than shrubs (Grass W - M. spinosum),
• Maximum LFMC in grasses is lower than shrubs (Grass W - S. bracteolactus)
• Maximum LFMC in grasses is lower than shrubs (Grass E - M. spinosum)
• Maximum LFMC in grasses is lower than shrubs (Grass E - S. bracteolactus)
• There are not difference between the two shrub species (M. spinosum - S. bracteolactus).

# 2.3.2 - Parameter "w" ----
contrast(emmeans(M1, ~leaf.type, param = "w"), "pairwise")

 contrast                       estimate   SE  df t.ratio p.value
 Grass W - Grass E                  20.7 2.59 226   7.981  <.0001
 Grass W - M. spinosum             -31.7 7.75 226  -4.085  0.0004
 Grass W - S. bracteolactus        -27.0 8.07 226  -3.342  0.0053
 Grass E - M. spinosum             -52.4 6.62 226  -7.914  <.0001
 Grass E - S. bracteolactus        -47.7 6.84 226  -6.974  <.0001
 M. spinosum - S. bracteolactus      4.7 6.72 226   0.699  0.8974

P value adjustment: tukey method for comparing a family of 4 estimates 

• Minimum LFMC in grasses from the W site is higher than grasses from the E site (Grass W - Grass E)
• Minimum LFMC in grasses is lower than shrubs (Grass W - M. spinosum)
• Minimum LFMC in grasses is lower than shrubs (Grass W - S. bracteolactus)
• Minimum LFMC in grasses is lower than shrubs (Grass E - M. spinosum)
• Minimum LFMC in grasses is lower than shrubs (Grass E - S. bracteolactus)
• There are not difference between the two shrub species (M. spinosum - S. bracteolactus).

3) NONLINEAR FIXED-EFFEC

Response function: lfmc = (A - w) / (1 + exp((m - time)/s))) + w

# 3.1.1 - Base nonlinear model ----
nl.fe.0 <- nls(lfmc ~ ((A - w) / (1 + exp((m - time) / s))) + w,
  start = c(A = 15, m = 30, s = -17, w = 5),
  data = lfmc
)
# here, time is the only predictor variable
nl.fe.0

Nonlinear regression model
  model: lfmc ~ ((A - w)/(1 + exp((m - time)/s))) + w
   data: lfmc
     A      m      s      w 
194.97  29.66 -21.17  27.93 
 residual sum-of-squares: 1029759

Number of iterations to convergence: 8 
Achieved convergence tolerance: 9.57e-07

b <- coef(nl.fe.0) # coefficients of the base model
b

        A         m         s         w 
194.97252  29.66004 -21.17131  27.92798 

# 3.1.2 - Leaf type fixed effect ----
nl.fe.1 <- nls(lfmc ~ ((A[leaf.type] - w[leaf.type]) / (1 + exp((m - time) / s))) + w[leaf.type],
  start = list(A = rep(b[1], 4), m = 25, s = -10, w = rep(b[4], 4)),
  data = lfmc
) # and fit is made using the base model's coefficients as start values
nl.fe.1

Nonlinear regression model
  model: lfmc ~ ((A[leaf.type] - w[leaf.type])/(1 + exp((m - time)/s))) +     w[leaf.type]
   data: lfmc
     A1      A2      A3      A4       m       s      w1      w2      w3      w4 
 56.726  92.483 298.835 318.207  30.431 -20.912  27.172   2.795  44.528  38.298 
 residual sum-of-squares: 67786

Number of iterations to convergence: 5 
Achieved convergence tolerance: 5.958e-06

b1 <- coef(nl.fe.1) # coefficients of the model with leaf type and time as predictor variables
b1

        A1         A2         A3         A4          m          s         w1         w2         w3         w4 
 56.725716  92.482865 298.834954 318.207067  30.430747 -20.911576  27.171667   2.795385  44.527958  38.297817 

# 3.1.3 - Plot fixed effect ----
nl.fe.2 <- nls(lfmc ~ ((A[group] - w[group]) / (1 + exp((m - time) / s))) + w[group],
  start = list(A = rep(c(b1[1], b1[2], b1[3], b1[4]), 3), m = 25, s = -10, w = rep(c(b1[7], b1[8], b1[9], b1[10]), 3)),
  data = lfmc
) # the model is fit using "b1" as the start values
nl.fe.2

Nonlinear regression model
  model: lfmc ~ ((A[group] - w[group])/(1 + exp((m - time)/s))) + w[group]
   data: lfmc
     A1      A2      A3      A4      A5      A6      A7      A8      A9     A10     A11     A12       m       s 
 53.683  54.282  61.906 111.562 314.131 333.831  77.465 298.161 321.433  87.553 279.437 292.443  30.413 -20.696 
     w1      w2      w3      w4      w5      w6      w7      w8      w9     w10     w11     w12 
 28.369  27.480  25.962   0.890  43.524  41.223   6.091  26.610  39.495   2.265  66.661  38.222 
 residual sum-of-squares: 53914

Number of iterations to convergence: 6 
Achieved convergence tolerance: 7.249e-06

IC_tab(nl.fe.0, nl.fe.1, nl.fe.2)

# Residual checking:
par(mfrow = c(1, 1))
plot(nl.fe.2) # heterocedasticity in the residuals

hist(resid(nl.fe.2), main = "", xlab = "Residuals") # slightly skew distribution

qqnorm(resid(nl.fe.2))
qqline(resid(nl.fe.2))

This model assumes:

• “A” and “w” vary with leaf type and plot (fixed-effects).
• “m” and “s” are unique for all the plots and leaf types.
• Variance is homogeneous
• There is no temporal dependence
• Residuals following a normal distribution

M2 <- nl.fe.2
summary(M2)

Formula: lfmc ~ ((A[group] - w[group])/(1 + exp((m - time)/s))) + w[group]

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
A1    53.683      8.602   6.240 2.20e-09 ***
A2    54.282      8.636   6.285 1.72e-09 ***
A3    61.906      8.881   6.971 3.59e-11 ***
A4   111.561     13.291   8.394 5.65e-15 ***
A5   314.131     24.840  12.646  < 2e-16 ***
A6   333.831     26.640  12.531  < 2e-16 ***
A7    77.465     10.997   7.044 2.34e-11 ***
A8   298.161     24.917  11.966  < 2e-16 ***
A9   321.433     25.765  12.475  < 2e-16 ***
A10   87.553     11.744   7.455 2.01e-12 ***
A11  279.437     20.829  13.416  < 2e-16 ***
A12  292.443     24.178  12.095  < 2e-16 ***
m     30.413      2.896  10.504  < 2e-16 ***
s    -20.695      3.606  -5.739 3.11e-08 ***
w1    28.369      6.537   4.340 2.17e-05 ***
w2    27.480      6.548   4.197 3.93e-05 ***
w3    25.962      6.633   3.914 0.000121 ***
w4     0.890      8.173   0.109 0.913389    
w5    43.524     13.567   3.208 0.001535 ** 
w6    41.223     14.429   2.857 0.004686 ** 
w7     6.091      7.262   0.839 0.402467    
w8    26.610     13.604   1.956 0.051725 .  
w9    39.495     14.009   2.819 0.005252 ** 
w10    2.265      7.551   0.300 0.764455    
w11   66.661     11.478   5.808 2.19e-08 ***
w12   38.222     13.031   2.933 0.003708 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 15.62 on 221 degrees of freedom

Number of iterations to convergence: 6 
Achieved convergence tolerance: 7.249e-06

4 - LINEAR MIXED-EFFECTS MODEL (M3)

Response function: lfmc = β0 + β1xTime + β2xGE + β3xSM + β4xSS β5xGExTime + β6xSMxTime + β7xSSxTime

where GE, SM, and SS are dummy variables (0 or 1) created to indicate differences between the base level of “leaf type”, i.e., grasses from the W site [GW], and the remaining levels.

4.1.1 Base linear mixed-effect mode

# 4.1.1 - Base linear mixed-effect model ----
l.me <- lme(lfmc ~ time * leaf.type, random = ~ 1 | plot, data = lfmc)

# Residual checking:
plot(l.me) # heterocedasticity in the residuals

plot(ACF(l.me, resType = "normalized", na.action = na.omit), alpha = 0.05, grid = TRUE) # temporal correlation is observed at lag 1

# 4.1.2 - Variance modelling ----
l.me.vm <- update(l.me, weights = varIdent(form = ~ 1 | leaf.type))

# Residual checking:
plot(l.me.vm) # variances there seem to be well modeled

AIC(l.me, l.me.vm) # modeling the variance improves the model fit

# 4.1.3 - Temporal correlation modelling ----
l.me.vm.tm <- update(l.me.vm, correlation = corARMA(p = 1, q = 1, form = ~1))
plot(ACF(l.me.vm.tm, resType = "normalized", na.action = na.omit), alpha = 0.05, grid = TRUE) # the temporal dependence is modeled

AIC(l.me.vm, l.me.vm.tm) # modeling the temporal correlation improves the model fit

plot(l.me.vm.tm)

plot(ACF(l.me.vm.tm, resType = "normalized", na.action = na.omit), alpha = 0.05, grid = TRUE)

hist(resid(l.me.vm.tm, type = "normalized"), main = "", xlab = "Residuals")

qqnorm(resid(l.me.vm.tm, type = "normalized"))
qqline(resid(l.me.vm.tm, type = "normalized"))

# 4.1.4 - Evaluation of the fixed-effects ----
l.me.ml <- update(l.me.vm.tm, method = "ML") # the model is fitted by Maximum likelihood
l.me.ml.1 <- update(l.me.ml, ~ . - time:leaf.type) # the model is reduced removing the interaction "time x leaf type"
IC_tab(l.me.ml, l.me.ml.1) # the AIC comparison suggests the interaction term to be important

M3 <- l.me.vm.tm
summary(M3)

Linear mixed-effects model fit by REML
  Data: lfmc 
       AIC     BIC    logLik
  1998.483 2050.63 -984.2416

Random effects:
 Formula: ~1 | plot
        (Intercept) Residual
StdDev:    3.750488 7.842661

Correlation Structure: ARMA(1,1)
 Formula: ~1 | plot 
 Parameter estimate(s):
      Phi1     Theta1 
 0.9076001 -0.7567060 
Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | leaf.type 
 Parameter estimates:
         Grass W      M. spinosum S. bracteolactus          Grass E 
        1.000000         2.963571         3.123740         1.070897 
Fixed effects:  lfmc ~ time * leaf.type 
                                   Value Std.Error  DF    t-value p-value
(Intercept)                     50.14877  3.478211 234  14.417977       0
time                            -0.27535  0.048114 234  -5.722819       0
leaf.typeGrass E                24.40614  4.908714 234   4.972002       0
leaf.typeM. spinosum           194.85292  8.539072 234  22.818980       0
leaf.typeS. bracteolactus      209.78398  8.838810 234  23.734415       0
time:leaf.typeGrass E           -0.56367  0.071704 234  -7.861042       0
time:leaf.typeM. spinosum       -2.02767  0.152372 234 -13.307359       0
time:leaf.typeS. bracteolactus  -2.29338  0.160937 234 -14.250193       0
 Correlation: 
                               (Intr) time   lf.tGE lf.M.s lf.S.b tm:.GE t:.M.s
time                           -0.557                                          
leaf.typeGrass E               -0.709  0.395                                   
leaf.typeM. spinosum           -0.407  0.227  0.654                            
leaf.typeS. bracteolactus      -0.394  0.219  0.653  0.666                     
time:leaf.typeGrass E           0.374 -0.671 -0.592 -0.430 -0.429              
time:leaf.typeM. spinosum       0.176 -0.316 -0.299 -0.742 -0.409  0.546       
time:leaf.typeS. bracteolactus  0.166 -0.299 -0.318 -0.441 -0.743  0.562  0.567

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-2.3449313 -0.6438200 -0.1011963  0.6158869  3.3469739 

Number of Observations: 247
Number of Groups: 6 

M4 <- lm(lfmc ~ time * group, data = lfmc)
summary(M4)

Call:
lm(formula = lfmc ~ time * group, data = lfmc)

Residuals:
    Min      1Q  Median      3Q     Max 
-49.673  -7.042   0.466   7.002  66.332 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)          48.79275    6.44309   7.573 9.60e-13 ***
time                 -0.24478    0.13305  -1.840  0.06714 .  
groupGW plot 5        0.49817    9.11190   0.055  0.95645    
groupGW plot 6        6.14984    9.11190   0.675  0.50042    
groupGE plot 1       40.32675    9.49655   4.246 3.19e-05 ***
groupSM plot 1      212.68278    9.11190  23.341  < 2e-16 ***
groupSS plot 1      227.42468    9.11190  24.959  < 2e-16 ***
groupGE plot 2       14.39850    9.49655   1.516  0.13089    
groupSM plot 2      195.54780    9.11190  21.461  < 2e-16 ***
groupSS plot 2      217.12857    9.11190  23.829  < 2e-16 ***
groupGE plot 3       21.68513    9.49655   2.283  0.02334 *  
groupSM plot 3      189.79613    9.49655  19.986  < 2e-16 ***
groupSS plot 3      192.53490    9.49655  20.274  < 2e-16 ***
time:groupGW plot 5  -0.01905    0.18816  -0.101  0.91944    
time:groupGW plot 6  -0.10231    0.18816  -0.544  0.58718    
time:groupGE plot 1  -0.80129    0.19357  -4.139 4.94e-05 ***
time:groupSM plot 1  -2.36256    0.18816 -12.556  < 2e-16 ***
time:groupSS plot 1  -2.55767    0.18816 -13.593  < 2e-16 ***
time:groupGE plot 2  -0.43459    0.19357  -2.245  0.02574 *  
time:groupSM plot 2  -2.34722    0.18816 -12.474  < 2e-16 ***
time:groupSS plot 2  -2.45552    0.18816 -13.050  < 2e-16 ***
time:groupGE plot 3  -0.56657    0.19357  -2.927  0.00378 ** 
time:groupSM plot 3  -1.82099    0.19357  -9.407  < 2e-16 ***
time:groupSS plot 3  -2.16847    0.19357 -11.202  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 16.51 on 223 degrees of freedom
Multiple R-squared:  0.9586,    Adjusted R-squared:  0.9544 
F-statistic: 224.7 on 23 and 223 DF,  p-value: < 2.2e-16

# Residual checking:
par(mfrow = c(2, 2))
plot(M4) # a clear pattern in the residuals is observed (model assumptions are not met)

# 6 - NULL MODEL (M5) ######################################################################################

M5 <- lm(lfmc ~ 1, data = lfmc)
summary(M5)

Call:
lm(formula = lfmc ~ 1, data = lfmc)

Residuals:
   Min     1Q Median     3Q    Max 
-88.75 -58.05 -31.45  52.66 231.06 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   95.645      4.919   19.45   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 77.3 on 246 degrees of freedom

IC_tab(M2, M3, M4, M5) # , weights = T, delta = TRUE, base=T, sort = TRUE) # the nonlinear mixed-effects model is clearly the best

IC_tab(M1)

par(mfrow = c(2, 2), cex = 0.65) # (Figure 5)
v1 <- c(1, 2, 3, 4)
v2 <- c("GW", "GE", "SM", "SS")

# Nonlinear mixed-effects model (M1)
rM1 <- resid(M1, type = "normalized")
fM1 <- fitted(M1)
plot(fM1, rM1, xlab = "Fitted LFMC (%)", ylab = "Residuals", ylim = c(-3, 3))
abline(a = 0, b = 0, lwd = 2)

boxplot(rM1 ~ lfmc$leaf.type, ylab = "Residuals", xlab = "Leaf type", xaxt = "n", ylim = c(-3, 3))
axis(side = 1, at = v1, labels = v2, las = 1)

plot(rM1 ~ lfmc$time, ylab = "Residuals", xlab = "Time (days)", ylim = c(-3, 3))
abline(a = 0, b = 0, lw = 2)

qqnorm(rM1, main = "", xlab = "Theoretical quantiles", ylab = "Sample quantiles")
qqline(rM1)

# Nonlinear fixed-effects model (M2)
par(mfrow = c(2, 2), cex = 0.65) # (Figure 5)

rM2 <- resid(M2)
fM2 <- fitted(M2)
plot(fM2, rM2, xlab = "Fitted LFMC (%)", ylab = "Residuals")
abline(a = 0, b = 0, lwd = 2)

boxplot(rM2 ~ lfmc$leaf.type, ylab = "Residuals", xlab = "Leaf type", xaxt = "n")
axis(side = 1, at = v1, labels = v2, las = 1)

plot(rM2 ~ lfmc$time, ylab = "Residuals", xlab = "Time (days)")
abline(a = 0, b = 0, lw = 2)

qqnorm(rM2, main = "", xlab = "Theoretical quantiles", ylab = "Sample quantiles")
qqline(rM2)

par(mfrow = c(2, 2), cex = 0.65)

# Linear mixed-effects model (M3)
rM3 <- resid(M3, type = "normalized")
fM3 <- fitted(M3)
plot(fM3, rM3, xlab = "Fitted LFMC (%)", ylab = "Residuals", ylim = c(-3, 3))
abline(a = 0, b = 0, lwd = 2)

boxplot(rM3 ~ lfmc$leaf.type, ylab = "Residuals", xlab = "Leaf type", xaxt = "n", ylim = c(-3, 3))
axis(side = 1, at = v1, labels = v2, las = 1)

plot(rM3 ~ lfmc$time, ylab = "Residuals", xlab = "Time (days)", ylim = c(-3, 3))
abline(a = 0, b = 0, lw = 2)

qqnorm(rM3, main = "", ylab = "Sample quantiles", xlab = "Theoretical quantiles")
qqline(rM3)

par(mfrow = c(2, 2), cex = 0.65)

# Clasical regression (M4)
rM4 <- resid(M4)
fM4 <- fitted(M4)
plot(fM4, rM4, xlab = "Fitted LFMC (%)", ylab = "Residuals")
abline(a = 0, b = 0, lwd = 2)

boxplot(rM4 ~ lfmc$leaf.type, ylab = "Residuals", xlab = "", xaxt = "n")
axis(side = 1, at = v1, labels = v2, las = 1)

plot(rM4 ~ lfmc$time, ylab = "Residuals", xlab = "Time (days)")
abline(a = 0, b = 0, lw = 2)

qqnorm(rM4, main = "", xlab = "Theoretical quantiles", ylab = "Sample quantiles")
qqline(rM4)

—-